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ABSTRACT 

Flutter boundaries, as well as flutter limit cycle amplitudes, fre- 
quencies and stresses were computed for a panel of length-width ratio 
4.48 exposed to applied in-plane and transverse loads. The Mach number 
range was 1.1 to 1.4 . The method used involved direct numerical 
integration of modal equations of motion derived from the nonlinear plate 
equations of von Karman, coupled with linearized potential flow aero- 
dynamic theory. The results obtained were compared to experimental data 
reported in Ref. 5. 

The flutter boundaries agreed reasonably well with experiment, except 
when the in-plane loading approached the buckling load. Structural damping 
had to be introduced, to produce frequencies comparable to the 
experimental values. Attempts to compute panel deflections or stress at a 
given point met with limited success. There is some evidence, however, 
that deflection and stress maxima can be estimated with somewhat greater 


accuracy . 
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NOMENCLATURE 


a 


= panel length 


a 


n 


* modal amplitude 


b 


= panel width 


B...„ = nonlinear elastic terms 

ljkfc 
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g s “ D 


(P_ha 4 
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1/2 


H.., I., 

n ji 


K = a) 


M 
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N 
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= speed of sound 
* panel bending stiffness 

= structural damping factor 
= dimensionless structural damping factor 
= panel thickness 

= aerodynamic admittance functions (Eq. 10) 

= dimensionless flutter frequency 
= Mach Number 

= number of modes used (Eq. 3) 

R x /rR x)buckle 
= pressure 


q - j PU 2 

V 

R x 

s = (X /y) 1 / 2 r 


= dynamic pressure 

= generalized aerodynamic force (eq. 10) 
= streamwise applied in-plane load 
= dimensionless aerodynamic time 
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= time 
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x,y 


Greek 


6. . 


Ap 


- ■ ¥ 


pU2j 
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4 = P a /P m h 
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a , o 
x’ y 


T = t 


' o r 


ip 


ha 4 
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r m 


CD 


= flow velocity 
~ panel deflection 
* coordinates in plane of plate 
= coordinate normal to plate 

= Kronecker delta 
= static pressure differential 
= dimensionless static pressure differential 
= dimensionless flow dynamic pressure 
= dimensionless flow density 
= Poisson's ratio 
= flow density 

= panel density 
= panel stresses 
= dimensionless time 
= velocity potential 
= Airy stress function 
= modal function (Eq. 3) 

= flutter frequency 
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L. INTRODUCTION 

It is now well established that panel flutter is not, in many cases, 
an immediately destructive vibration* Hence flutter may be tolerated if it 
can be established that the flutter amplitude is sufficiently small and the 
duration of flutter sufficiently short. Unfortunately linear structural 
and aerodynamic theory is incapable of determining flutter amplitudes. Only 
by including the important panel nonlinearities can the flutter amplitude 
be established. Recently, methods have been developed at Princeton for 
analyzing the large amplitude oscillations of a fluttering plate. In the 
investigation reported here, these methods were used to calculate the 
flutter behavior of a panel exposed to a static pressure differential (that 
is, an applied transverse pressure load), and to applied in-plane com- 
pressive loads comparable to the buckling load of the panel. 

The panel length-width ratio (4.48), and the Tange of flow Mach number 
(1.1 to 1.4) were selected to allow comparison with the results of wind 
tunnel tests reported in Reference S. These tests were in turn motivated by 
a desire to investigate the flutter behavior of certain panels mounted on 

5 

the forward shirt of the S IV-B stage of the Saturn V launch vehicle. During 
these tests, the frequency and amplitude of the panel motion (if any) were 
measured as the tunnel dynamic pressure was increased. The tests were carried 
out at various values of test section Mach number, panel static pressure dif- 
ferential, and applied in-plane load. By this method both the panel flutter 
boundaries (lowest dynamic pressure at which flutter occurrred) and the severity 
of the post flutter motion were determined. 

The calculations described herein were carried out for the same range of 
parameters as used in Ref. S. The method used involves the direct numerical 
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integration of a set of nonlinear differential equations for the panel 
motion, derived from an approximate modal solution of the von Karman 
nonlinear plate equations. Because of the range of Mach numbers involved, 
the popular quasi-steady or piston theory expressions for the aerodynamic 
pressure on the panel were not applicable. Instead the full linearized 
inviscid, potential flow theory was employed. 

So far as is known, the work reported herein constitutes the first 
attempt at predicting theoretically the severity of flutter of a low 
aspect ratio stressed panel in the critical low supersonic Mach number 
range. 



- 3 - 


II. THEORETICAL DEVELOPMENT 


The equations of motion for a three dimensional plate, von Kartnan's 
large deflection equations are 


DV“w = 92 * * 2w 
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( 2 ) 


where w is the plate deflection and $ is the Airy stress function. is 

a structural damping parameter. The reason for including structural damping 
will be discussed later. Equation (2) and the first three terms on the right 
hand side of equation (1) constitute the nonlinear elastic coupling between 
out-of-plane bending and in-plane stretching that ultimately limits the 
amplitude of flutter. 

Equations (1) and (2) are reduced to a set of simultaneous nonlinear 
differential equations by Galerkin’s method. The transverse displacement w 
is expressed as a linear combination of modal functions that satisfy the 
appropriate boundary conditions at the edge of the plate (in this case, those 
for a clamped plate) : 

N 

w/h = i a (t) ^ (x/a) ip (y/b) 

tn— 1 1 


cos (m - 1)? - cos (m + 1)? 
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As is described in greater detail in Refs. 1-4, $ is determined by solving 

equation (2) with expression (3) inserted for w . The boundary conditions 
satisfied by $ on the plate edges (the so-called in-plane boundary condi- 
tions) depend on the design of the panel support structure. In Refs. 3 and 4 
methods of handling situations corresponding to either complete restraint 
(no in-plane motion permitted at the edges) or zero restraint (in-plane 
stresses zero at the edges) are discussed. It is not generally feasible 
to distinguish between these two alternate sets of boundary conditions before- 
hand by analyzing the panel support structure (and in fact most practical 
structures would create a degree of restraint somewhere between the two 
extremes), so both sets are retained in the developments that follow. With 
$ determined, equation (1) is satisfied in the Galerkin sense by computing 
the integral average of equation (1) weighted successively by each of the 
modal functions ip . (x/a) ^^(y/b) in expression (3) and setting the result 
to zero. The resulting system of equations is (in nondimens ional form): 


E 

j 


S. . (a. + g a.) + E C. .a. 
3 B s y j 30 j 


+ Z II 
j k A 


ijkA a j\ a £ + 


* 

X E 



AP 6 «* 0 

“ H 


(4) 


The matrices S and C are the familiar modal mass and elastic stiffness 
matrices of linear vibration theory. C contains the applied streamwise 
in-plane tension as a parameter; when R^ decreases below a critical 

negative value, the plate buckles. The fourth order array B contains the 
nonlinear terms corresponding to the coupling between in-plane stretching and 
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and out-of-plane bending referred to previously. Explicit expressions for 

all of these terms are contained in Ref. 4. 

The generalized aerodynamic forces Q. . are defined as 

J i 




V '■ 


0 0 


Pj ~Pa 


PIP 2 


^(x/aH (y/b) 


(5) 


where p^ is the pressure on the plate caused by an arbitrary deflection in 
the jth mode: 


w = aj (t) ipj (x/a)^ (y/b) 


pj is given by 




z=0 


where the velocity potential 4 1 must satisfy 
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subject to the boundary conditions 
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off plate 


(9) 
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The boundary value problem defined by Eqs . (6-9) has been solved in Ref. 7, 
where it is shown that 
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da. 


v (a.D.. + -r± S..) 

M j ji ds Ji' 


is 

| a., (a) H_..(a) d0 
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da (a) 

d^ I^Cs-^da 


( 10 ) 


with 

s = (\*/y) l/2 T (11) 

See Ref. 7 or Appendix B of Ref. 4 for evaluations of D.., S.., H..(s), 

ji’ ji’ ji w ’ 

and I„(s). (Beware of slight notational differences between the two.) 

These functions depend parametrically on M and a/b , but not explicitly 
* 

on X and y . If the integrals in (10) are deleted, the are those 

given by "piston theory"; that is by a direct substitution of the well known 
expression 

pU 2 f 3w 1 3w 

P " P « = M [ax + U 3t 

into equation (1) . 

Equations (4) and (10) are combined to form a set of coupled nonlinear 

ordinary integral-differential equations in time, t . The solution 

* 

procedure is to specify X , y, M, Ap, a/b, R , g , and to determine the 

X s 

modal amplitudes by numerical integration. Given the a^fr), the deflection 
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w/h or stresses a , o at any selected point on the panel may be cal- 

x y 

culated in a straightforward manner. The computer programs used to 

carry out these various procedures are listed in the Appendix. 

These routines are modified and improved version of the programs listed in 


Ref. 4 ♦ 
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III, COMPUTATIONAL CONSIDERATIONS 

The considerations of this section relate to the manner in which the 
computations were arranged and carried out, and to the way in which the 
results obtained have been displayed. They have been dictated both by the 
nature of the wind tunnel experiments reported in Ref. 5, and by the 
necessity of using a very large number of modes (12 in most cases) in order 
to properly represent the behavior of the low aspect ratio panel being 
studied. 

In order to save computer time (and hence expense) it was found 
useful to divide the computations into four distinct steps. These are: 

1) Computation of the nonlinear terms. Only the plate length width 
ration a/b, the Poisson’s ratio v , and the in-plane boundary conditions 
need be specified in order to determine B. Since the results were to be 
compared with the data of Ref. 5, only one value of a/b (= 4.48) and 

v (* 0.3) were employed. Hence only two sets of nonlinear terms, cor- 
responding to complete and zero in-plane edge restraint, were required. 

These were computed at the outset, and stored on magnetic tape. 

2) Computation of the aerodynamic admittance functions H^(a/b, M, s) 

H„(a/b, M, s) (see Eq.' 10). As indicated, these quantities depend on the 
panel length-width ratio and the flow Mach number as well as the dimen- 
sionless aerodynamic time s. Since only four values of M were studied, 
it was found worthwhile to compute and I ^ j beforehand as well 

(distinct sets of values for each of M * 1.1, 1.2, 1.3, and 1.4). They 
also were stored on magnetic tape. 

3) Numerical integration of the panel equations of motion. This 
operation uses as inputs the data stored from steps 1) and 2) above. 
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as well as specification of X , P» M, Ap, R x> and g g . Interest 

centers on the amplitude and frequency of the flutter limit cycle at 
a given point on the panel: 

w/h) - f(x\ y, M, Ap, R x , g s ) 

K - g( " ) 

(The cross-stream in-plane load R was zero in the experiments of 

Ref. 5. and so was assigned the same value in the present study.) 

* 

The dimensionless flow dynamic pressure X and flow density y 
are related through the flow velocity: 

X 
P 


p ha 2 
m 


D 


U 2 


The quantity in brackets is uniquely defined by the geometric and material 
specifications of the panel being studied. Furthermore, in a continuous 
flow wind tunnel the flow velocity U is determined by the test section 
Mach number M, and the stagnation temperature T q in the tunnel 
settling chamber: 

u = C RT o ) Kl 

''O'* 

M 2 T _ M 2 

T o 1+3— M 2 


The stagnation temperature is held constant during tunnel operation, 
so U is determined solely by the test section Mach number. 
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It is therefore convenient to display the results of the flutter 
amplitude and frequency as functions of q, the (dimensional) dynamic 
pressure rather then as functions of both X and u independently: 

w/h) = F(q, M, Ap, N x , g g ) 

K = G(q, M, Ap, N x , g g ) 

Since the non-dimensionalization of R x is arbitrary, it has been 
replaced here by the ratio of R x to its buckling value: 


N 


x 


^x^x^buckle 


By extrapolating to w/h) 0, it is possible to determine the critical 
or flutter dynamic pressure q £ and the flutter frequency K £ : 

q £ = F^ (M, Ap, N x , g s ) 

K £ = G^M, p, N x , g s ) 

4) Panel stresses during flutter. In the theory of thin plates, 
normal stresses vary linearly across the plate thickness. The extreme 
values of stress occur on the upper and lower surfaces of the panel, e.g. 

x x^ms x'b 

where the + and - sign apply to the upper and lower surfaces, 
respectively. A similar equation holds for a . The bending stress 
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a is proportional to the local curvature of the plate, and is obtained 
from the modal amplitudes a n by differentiating Eq. (3) for w/h . On 
the other hand, the middle surface or in-olane stress a ) obtained by 
differentiating the Airy stress function $ of Eq. (2). As such the 
in-plane stresses depend not only on the plate deflection w(x/a, y/b,x) 
but also on the in-plane boundary conditions satisfied at the edges of 
the plate. The computer program listed in the Appendix uses the modal 
amplitudes a^t) from step 3) to calculate the in-plane or middle- 
surface stress for a panel completely restrained at its edges . Since 
not many flutter calculations were made for the zero edge restraint 
case, an equivalent program for computing the middle surface stresses 
in such panels was not written. 
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IV. NUMERICAL RESULTS 


Free Panel Vibrations 

In order to explore the extent to which the theoretical model 
employed mirrors the elastic behavior of the panel, independently 
of the flutter results, panel natural frequencies were computed as 
a function of applied static transverse and in-plane loading. The 
transverse load was equivalent to a pressure differential between 
the two faces of the panel . 

The computations were carried out by integrating the modal equations 
* 

(4) (with X * p =0, g g > 0) to determine the equilibrium panel 
deflection under the assumed loading, and then linearizing the equations 
about that deflection. The natural frequencies were determined numeri- 
cally from these linearized equations by solving a classical eigenvalue 
problem. Representative results are shown in Figures 1 through 4, 
along with comparable experimental data from Ref. 5. 

Figures 1 through 3 show the behavior of modes 1 , 2 , and 6 under 
a transverse pressure loading. In each case calculations were made 
assuming both zero and complete in-plane edge restraint. (The edges 
of a plate with zero in-plane restraint are free to move in the plane 
of the plate in response to transverse plate motions, while the edges 
of a plate with complete in-plane restraint are prevented from making 
any such movement .) In all three figures there is a systematic dis- 
crepancy at zero pressure load. Part of this difference is attributable 
to imperfect convergence of the solution, but probably not all. The 
calculated frequencies included in Ref. 5 (Table II) show a similar deviation 
from the experimental results. Of greater interest, however, are the amounts 
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by which the various frequencies increase when a pressure load is applied. 
For the lower modes, the assumption of complete edge restraint provides 
the best agreement with experiment, while for the higher modes, zero 
edge restraint works best. 

Figure 4 shows calculated and experimental results for the behavior 
of the ninth mode under a compressive in -plane load. Both calculated 
frequencies drop off much more near the buckling load than does the 
experimental curve. This may be due to the presence of imperfections 
in the plate, such as a slight initial curvature or waviness. 

It was not possible, on the basis of these results, to eliminate one 
in-plane boundary condition from further consideration. Therefore flutter 
calculations were carried out for both cases, although a shortage of time 
and money limited the number of zero edge restraint runs that could be 
carried out. 

Flutter Calculations - General Nature of Solutions Obtained 

The flutter limit cycle was determined by integrating the modal 
equations (4) until a periodic motion was found. Experience indicates 
that the initial conditions used to start the integration do not affect 
the amplitude or frequency of the limit cycle, at least for N < 1 . 
Because of the large length-width ratio (a/b = 4.48) employed, at least 
twelve modes were required to obtain an acceptable degree of convergence. 
Furthermore, the transient portion of the solution survived for the 
equivalent of many cycles of the ultimate limit cycle motion. (Neither 
of these statements apply for smaller values of a/b.) As a result, the 
numerical integration turned out to be costly in terms of both computer 
memory storage area and computation time . 
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Initially all calculations were made without introducing any damping 
other than the aerodynamic damping implicit in the potential flow expression 
(10) for the generalized forces Qj^* However, it was found that for larger 
values of the dynamic pressure q the flutter frequency became very high 
(X 900 Hz), with the panel deflection being such that the 9th or 10th mode 
had the largest amplitude. Cunningham 8 has shown that the flutter frequency 
and mode shape are both critically sensitive to the amount of structural 
damping present. Therefore, structural damping was introduced into the 
panel equations of motion (4) in order to suppress the high frequency 
flutter. 

The structural damping present in the actual panels used in Ref. 5 
has not been measured to date. Moreover, if the lower flutter frequency 
found experimentally is indeed due to the presence of additional damping, 
the source of that damping need not be structural. It may well be caused by 
the boundary layer in the airflow over the panel. The capability for dealing 

9 

with boundary layer effects does exist, but at least for the present 
the technique is not practical for flutter calculations involving the use 
of many structural modes . Hence the introduction of structural damping 
must be viewed as an essentially ad hoc procedure designed to eliminate 
a physically unrealistic aspect of the flutter behavior. 

From a mathematical standpoint, there are many forms of structural 
damping that can be introduced into the plate equations to describe 
non-elastic behavior. Of these the traditional and most popular choice 
is the (1 + ig) type, which is meaningless for non-sinuousoidal motion 
and is therefore unsuitable for nonlinear plate equations. The most 
common of the many types that can be used have the general form 
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3w 

at 


with n = 0, 1, 2, .... 

They differ, in the modal formulation used here, in the relative damping 
ratios given the various modes. Roughly speaking, the damping ratios 
increase as the mode number raised to the (n-l)st power. For the present 
work, n = 2 was selected because G s fits easily into the modal 

equations (4) , and because it provides greater damping in the higher 
modes whose motion it is intended to suppress. 

Most of the results that follow have been calculated with 
S s = .0001. This provides a damping ratio for the first 
mode of .025 (2.5% of critical). Cost limitations have made it impossible 
to present a systematic study of the influence of structural damping over 
the complete range of Mach number, static pressure differential, and applied 
in-plane load considered here. 

Flutter Boundaries 

Figures 5 and 6 show flutter boundaries as a function of Mach number 
for an unloaded panel (Ap = = 0) . Curves are shown for g s = 0. and 

g s * .0001. On the same figures are shown experimental data from Ref. 5. 

The data were obtained from several different panels (of indentical specifi- 
cation), and for two different boundary layer thicknesses, the thicker one 
being induced by inserting spring pins in the tunnel wall ahead of the panel, 
which was mounted flush to the wall. This caused the boundary layer thickness 
(as measured near the trailing edge of the panel) to increase roughly 7 to 
30% , depending on the tunnel Mach number and dynamic pressure.^ 
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Both the theoretical and experimental results show a gradual decrease 
in flutter dynamic pressure with increasing Mach number, but the theory 
shows no minimum at M = 1.3 . The theoretical flutter boundary for 
g s = 0 . agrees best with the experimental results for the smooth wall 
boundary layer (Figure 5) while the damped flutter boundary agrees best 
with the results for the rough wall boundary layer (Figure 6) . This is 
the correct qualitative behavior, since the boundary layer introduces a 
damping effect that increases with boundary layer thickness. The quanti- 
tative agreement in Figure 6 is of course fortuitous, since the amount of 
structural damping introduced is arbitrary, and in any event the damping 
is of structural origin in the theory and aerodynamic in the experiment. 

Flutter frequencies are shown in Fig. 7. As mentioned previously, 
the frequencies for zero damping are unrealistically high, whereas those 
for g s = .0001 are comparable to the experimental results. Neither 
theory nor experiment shows much variation with Mach number. 

Figures 8 and 9 show calculated and experimental flutter boundaries 
for plates exposed to compressive in-plane loads. In both figures the 
qualitative behavior with N x is correct, although the rate of decrease 
in the flutter dynamic pressure is more rapid according to the theory. 

Near buckling (N^ = 1.0), the theoretical result becomes overly conserva- 
tive. The calculations of flutter boundaries near the buckling load is 
a difficult matter, since the plate behavior is then especially sensitive 
to the presence of small initial structural imperfections , to the damping 
effect of the boundary layer, and so on. The prediction of panel natural 
frequencies under in -plane loading suffers the same difficulty, as can 
be recalled from Fig. 4. 
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Figure 10 shows a limited set of calculations for a panel exposed 

to a static pressure differential. Both sets of in-plane boundary 

conditions (zero and complete edge restraint) are included. The line 

labeled ''Exp." is a derived curve taken from Fig. 43 of Ref. 5. The 

result for zero in-plane restraint shows the better agreement with 

experiment, in spite of the fact that the panels referred to in Ref. 5 

were carefully mounted in a massive supporting structure. This result 

is consistent with similar comparisons made previously involving panels 

3 10 

of smaller length-width ratio at higher Mach numbers. ’ 

Figs. 11, 12, and 13 contain flutter boundaries for panels exposed 

to combined loading (both Ap t 0. and N 0.). The theoretical results 

X 

are all for the case of complete edge restraint, and reflect the same 
behavior as exhibited in Fig. 10, namely, a lack of sensitivity to static 
pressure differential. Note, however, that the agreement between the 
slopes of corresponding pairs of flutter boundaries in Fig. 11 improves 
as the in-plane loading increases. 

It would be desirable to carry out calculations equivalent to those 
shown in Figs. 11-13 for the zero in-plane restraint case. 

Panel Displacement in Flutter 

A record of panel centerline deflection during approximately one period 
of the flutter oscillation is shown in Fig. 14. Structural damping (g g = .0001) 
was assumed in making the calculation; with no damping, many more zero cros- 
sings appear than are shown. The motion portrayed in Fig. 14 is quali- 
tatively similar to that reported in Ref. 5. (See especially Fig. 57 of 
that report). In particular, the panel deflection is largest near the 
trailing edge, but not markedly so, and the streamwise variation of the 
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deflection is elaborate, but with relatively few zero crossings at any 
given instant. The motion has a quasi wave-like character, since the 
zero crossings (points of zero deflection) move with time, and even 
appear and disappear. 

Relatively little panel displacement data was published in Ref. 5, 
and what was presented was limited to a case wherein the panel was 
buckled by the applied in-plane load. This situation is both the 
most important physically (since at a given dynamic pressure the panel 
deflection is maximized), and the most difficult to handle analytically. 

As mentioned previously, buckled panels are especially sensitive to effects 
that normally are either ignored entirely (such as initial imperfections), 
or handled very crudely (structural damping) . 

Panel displacements at three different streamwise locations (but 
2.5 inches off the centerline) are shown in Fig. 15. The 
streamwise locations of probes A, C, and F are shown in Fig. 14. At 
all three locations, the calculated displacements are considerably larger 
than their experimental counterparts. 

Stresses 

The bending stresses are generally considerably larger than the in-plane 
or axial stresses during flutter. Since the bending stresses are pro- 
portional to local panel curvature, the bending stress distribution generally 
resembles the panel deflection (see Fig. 16). Attempts to compute stresses 
at a given point on the panel are therefore hampered by the same difficulty 
encountered in calculating deflections: a small change in the flutter mode 

shape causes large errors in the stresses computed at that point. 
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Fig. 17 shows a comparison of calculated and experimental stresses 
as a function of flow dynamic pressure for a buckled panel. The open 
circles are the peak-to-peak bending stress in panel #6 at the location 
of gauge Bl, just off the center line of the panel near its trailing 
edge. The small triangles connected by straight lines are theoretical 
peak-to-peak stresses calculated for the same point, as well as for a 
point on the panel center-line, three quarters of the way back behind 
the leading edge. This latter location is the point where the maximum 
stress occurred, according to the theory. It should be noted that the 
applied in-plane load assumed for the calculations was only 73.5% of the 
theoretical buckling load, whereas the experiment was carried out with 
an in-plane load equal to 96% of the buckling load applied. As can be 
seen, the bending stress computed at the 3/4 chord point agrees better 
with the experimental result than does the value computed at the position 
of gauge Bl. If the stress measured at Bl is in fact the maximum stress 
that occurred, then the maximum stress is computed with greater accuracy 
than is the stress at Bl. 

Figs. 18 through 21 (each of which is divided into two parts ) show 
similar comparisons. In each figure part (a) shows theoretical and experi- 
mental bending stresses (Figs. 18 and 19) or axial stresses (Figs. 20 and 21) at 
the point referred to above. In part (b) the calculated data is the theoretical 
maximum stress on the panel, displayed alongside the same experimental 
data as in part (a). In general, the maximum stress (which is less 
sensitive to changes in the flutter mode shape) best reflects the experi- 
mental trends, at least for small N x . Near the buckling load, excessively 
large maximum stresses are predicted, presumably for the same reasons 


mentioned earlier. 



- 20 


V. CONCLUSIONS 

Flutter boundaries were computed numerically as functions of Mach 
number, in-plane loading, and static pressure differential. Comparison 
with experimental data indicate reasonably good correlation for Mach 
number and in-plane loading, except near the buckling load. The influence 
of static pressure differential depends on the in-plane boundary conditions 
assumed. Assuming zero restraint (edges free to move in plane) provided 
the best correlation with experiment, although not enough calculations were 
made to firmly establish this point. 

The flutter mode shapes calculated were in good qualitative agreement 
with experiment. The flutter frequency, however, proved to be sensitive 
to the amount of structural damping assumed. With no damping, the coupled 
flutter frequency was several times higher than the experimental value. 
Because flutter frequency is an important factor in determining panel 
fatigue life, future experimental programs should include a determination 
of panel damping. In addition, the theoretical model employed should 
be improved to include the damping effect of the boundary layer. 

Attempts to compute panel deflection and stresses during flutter 
met with limited success, particularly for buckled panels. There is 
some indication, however, that maximum deflections and stresses can be 
calculated with greater accuracy than deflections on stresses at a 
specific point. From a practical standpoint, knowledge of the maximum 
is sufficient to determine panel fatigue life; the stress distribution 
and mode shape are of lesser significance. 

Most of the difficulties encountered in this investigation stem 
from the large length-width ratio of the panel and the presence of 
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large in-plane loads. In this regard a wind tunnel test program using 
a carefully constructed high aspect ratio (a/b < 1) panel would be 
very helpful. Stream-wise buckling loads might well be included in 
the test program, but an extensive set of data should also be collected 
with little or no in-plane loading present. Such data would be of great 
help in assessing current theoretical methods without the perplexing 
but not fundamental difficulties associated with low aspect ratio and 
panel buckling. 
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APPENDIX 


Listing of Computer Programs 



LEVEL 21.6 ( MAY 72 ) 


OS/360 FORTRAN H 




COMPILER 

OPTIONS - N A 1 E= MAIN,OPT=02,LIN£CNT=58,SIZE=0000K, 




SOURCE, EBCDIC, NO LIST, NODECK , LOAD , NAP, NOEDI T, I D, NO XREF 



C 

PROGRAM TO COMPUTE NONLINEAR TERMS (COUPLING BETWEEN 



C 

IN-PLANE STRETCHING AND OUT OF PLANE BENDING) FOR CLAMPED 



C 

PLATE WITH COMPLETE IN-PLANE EDGE RESTRAINT 



C 

AB = PLAT3 L.ENGTH/NIDTH RATIO, NO=POISSON ' S RA.TIO 



C 

NV = # OF MODES 

ISN 

0002 


REAL NU 

ISN 

0003 


REAL II 

ISM 

0004 


DIMENSION 8(12,12,12, 12) 

IS N 

0005 


DNA(K) - (FLOAT (K) **2+16. *AB2) **2 

ISM 

0006 


DNB(K) = (FLOAT (K) **2+4. *AE2) **2 

ISM 

0007 


CSS (K, L, M) = .5* (CC (K,L-M) -CC (K,L+M) ) 

ISN 

0008 


CCC (K,L, M) = . i* (CC(K,L-M) + CC(K,L + N)) 

ISN 

0009 


GG (K, L , M) = CCC (K- 1 , L- 1 , M) -CCC (K - 1 ,L+ 1 , N) -CCC (K+ 1 , L- 1 , M) 




1 *CCC(K+1,L+1,M) 

ISN 

0010 


HH (K, L,M) = -PI >*FLOA T (L-1) **2* (CCC (K- 1 ,L- 1,M) 




1 -CCC (K + 1,L-1 ,rt) ) 




2 + Pl2*FLOAT (L + 1) **2* (CCC (K-1,L + 1,M) -CCC (K+1,L+1,N) ) 

ISN 

0011 


II(K,L,M) - -PI*FLOAT (L-1) * (CSS (K- 1 , L- 1 , M) -CSS (K+1,L-1,H) ) 




1 + P I*FLO AT (L + 1) * (CSS (K-1,L+1,M) -CSS(K+1,L+1,M) ) 

ISN 

0012 

1 

. FORMAT (130) 

ISN 

0013 

2 

FORMAT (131) 

ISN 

0014 


PI * 3. 14 159 

ISN 

0015 


PI2 = PI*PI 

ISN 

0016 


PI3 = PI2+PI 

ISN 

0017 


PI 4 = PI3*PI 

ISN 

0018 


READ (5,110) A3 , N U 

ISN 

0019 

110 

FORMAT (4E10. 3) 

ISN 

0020 


WRITE (6,1101) A3, NU 

ISN 

0021 

1 101 

FORMAT (2£ 12. 3) 

ISN 

0022 


WRITE (6,1) 

ISN 

0023 


READ (5, 123) NV 

ISN 

0024 


WRITE (6 , 120) NV 

ISN 

0025 

120 

FORMAT (15) 

ISN 

0026 


WRITE (6,1) 

ISN 

0027 


AB2 - A3* *2 

ISN 

0028 


AB4 = A 0**4 

ISN 

0029 


X = 1 . 0E- 12 

ISN 

0030 


DO 4 2 a = 1, NV 

ISN 

0031 


DO 42 N = 1, NV 

ISN 

0032 


MA = M-N 

ISN 

0033 


MB = M-N-2 

ISN 

00 34 


MC = M-N+2 

ISN 

0035 


MD = M + N. 

ISN 

0036 


ME = M + N-2 

ISN 

0037 


MP = M+N* 2 

ISN 

0038 


rm = n 

ISN 

0039 


RN = N 

ISN 

0040 


R MA = MA 

ISN 

0041 


RM B = MB 

ISN 

0042 


RMC = MC 

ISN 

0043 


RM D = MD 

ISN 

0044 


RM E = MS 

ISN 

0045 


RM F = MP 

ISN 

0046 


BA » -2. * (HM*RM0 + 2.) /i)NA (MA) 


- 1 -ft 



ISN 

0047 


BB = 

(RM— 1 .) *RNJ/DNA (MB) 

isn 

0048 


BC = 

(RM+ 1.) *RMO/DNA <MC) 

ISN 

0049 


BD = 

2.* (8M*R8A4-2. ) /UN A (MD) 

ISN 

0050 


BE = 

- (RM-1. ) *RMA/DNA (ME) 

ISN 

0051 


BF = 

- (RM4- 1. ) *R.1 A/DN A { MF) 

ISN 

0052 


BG = 

4.* (RM**2+1.)/DN3 (MA) 

ISN 

0053 


BH = 

-2.* (RM-1.) **2/DNB (MB) 

ISN 

0054 


BK = 

-2. * (RiH-1.) * * 2/D N B (MC) 

ISN 

0055 


BL = 

-4.*{RM**2*1.)/DNB(MD) 

ISN 

0056 


BM = 

2.*(RM-1.) **2/UNi3 (ME) 

ISN 

0057 


BN « 

2 . * (RH + 1 * ) ** 2/l)N B (MF) 

ISN 

0058 


BP = 

-2 . * RM/ (R1A**3+X) 

ISN 

0059 


B Q = 

(Rtl-1.) / {8.13**34- X) 

ISN 

0060 


BR = 

(RM*1.)/(R1C**3+X) 

ISN 

0061 


BS = 

2. *RM/R P. D* * 3 

ISN 

0C62 


BT = 

- (RM-1. )/(RM£**3*X) 

ISN 

0063 


BU = 

- (RM+ 1 . ) / 3 M F * * 3 

ISN 

0064 


DO 42 

I - 1 , N V 

ISN 

0065 


DO 42 

K = I, NV 

ISN 

0066 


IODD 

= I4-K + M4-N 

ISN 

0067 


IF (MOD(IODD,2) .NE.O) GO TO 38 

ISN 

0069 


BAA = 

(BA-B G) (I , K , M A) ♦ (BB-3H) *KH (I , K, MB) + 



1 

( BC- BK) *Hil (I , K, MC) * (8D-BL) *HH (I , K,MD) ♦ 



2 

(BE -BN) *hH (I, X, ME) ♦ (BF-BN) *HH (I, K ,MF) 

ISN 

0070 


BAA a 

-43.*(1.-NU**2) * AB2*PI2* E A A 

ISN 

0071 


B BB = 

RMA* ( E A - 3 i) *II(I,K,MA) + RMB* (BB-BH) *11 (I,K 



1 

♦ RMC* (BC-BK) *11 (I,K, 1C) + RtID* (BD-BL) *11 (I # K,f1D) 



2 

RME* (BE-Bd) *11(1, K,dE) + R MF* (BF- BN) *11 (I , K # MF) 

ISN 

0072 


BB B = 

-24. + ( 1.-N J**2) *AB2*PI 3*B8B 

ISN 

0073 


BCC = 

RHA**2* (BA-2.*B3*2.*BP) *GG (I,K,BA) 


1 ♦ RMB**2* (B3-2. *BH+2. *3Q) *GG (I,K,MB) 

3 ♦ RMC**2* (BC-2.+BK+2. *BR) *GG (I,K,MC) 

4 «• RMD**2* (BO-2. *3L*2. *BS) *GG (I, K, «D) 

5 + RHE**2*(Bi:-2. *384-2. *3T)*GG(I,K, ME) 

6 + RMF**2* (BF— 2. *BN+2 • * 3U) *GG (I , K, 8F) 


ESN 0074 


BCC = 12.* (1.-NU**2) * A 32*PI4*BCC 

ISN 0075 


B(I,K,8,N) = -A32* (BAA-2. *E EB+OCC) 

ISN 0076 


GO TO 42 

ISN 0077 

38 

B(I,K,M,N) = 0. 

ISN 0078 

42 

B (K,I,M,N) = B(I, K,M,N) 


C 

WRITE NONLINEAR TERMS ONTO TAPE 

ISN 0079 


WRITE (13) 3 

ISN 0080 


END FILE 10 

ISN 0081 


DO 45 I = 1,NV 

ISN 0082 


DO 45 J = 1 , N V 

ISN 0083 


DO 4 5 K = 1 , N V 

ISN 0084 

45 

WRITE (0,47) (3(I,J,K,L), L = 1,NV) 

ISN 0085 

47 

FORMAT (10212. 3) 

ISN 0086 


STOP 

ISN 0087 


END 


2 



LEVEL 21.6 ( KAY 72 ) 


OS/360 FORTRAN H 


ISN 0002 
ISN 0003 
ISN 0004 
ISN 0006 
ISN 0008 
ISN 0009 


COMPILER OPTIONS - MAME= MAIN, O?r=02, LINECNT=58 ,SIZE=0000K, 

SOURCE, EJCDIC, NOLIST, NODECK, LOAD, MAP, NOEDIT, ID, NOXBEF 
FU NCTION CC ( K , M) 

CC = 0. 

IF (K. EQ..M) CC = CC*. 5 
IF(K.EQ.-M) CC * CC*.'5 
RET0RN 
END 


3 



FORTRAN IV G LEVEL 21 


MAIN 


DATE = 73248 


0001 

0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 
001 1 
0012 

0013 

0014 


0015 


0016 

0017 

0018 

0019 

0020 
002 1 
0022 

0023 

0024 

0025 

0026 

0027 

0028 

0029 

0030 

0031 

0032 

0033 

0034 

0035 

0036 

0037 

0038 

0039 


C PROGRAM TO CJMPUT3 NONLINEAR TERMS (COUPLING 

C BETWEEN IN-PLANE STRETCHING AND OUT OF PLANE BENDING) 

C FOR PLATES WITH ZERO IN-PLANE EDGE RESTRAINT 

C AB = PL AT S LENGTH/W I DTH RATIO, NU=POISSON‘ S RATIO, 

C NV=* OF MODES 

REAL NU 

DIMENSION 3 (12,12, 12, 12) 

DIMENSION AA(3Q,3,30,J) ,V (14,30,3,14) ,G(30,3,14,14) 

DIMENSION AM (90,9 0) , AM W (90, 90) 

DIMENSION A VW (180) 

REAL*R AM, Art W 
REAL*8 AVW 

CSS (K , L , M) = . 5*(CC(K,L-M) -CC (K,L + M) ) 

CCC (K , L , M) = .5* (CC(K,L-M) ♦ CC(K,L+M)) 

F (I , J) = CC (I,J) 

FF (N, L) = F (N- 1 , L- 1) - P(N-1,L+1) - F(N+1,L-1) ♦ F(N + 1,L*1) 

FFF (K, M) = -FLOAT (M-1) **2* (F <M-1 ,K-1) - F (M- 1 , K+ 1 ) ) 

1 +• FLOAT (M+ 1) * * 2* (/ (rt+1,K-1) - F(M+1,K+1)) 

FFFF (K , M) * FLOAT (M-1) **4* (F (M-1 ,K-1) - F(M-1,K+1)) - 

1 FLOAT (M + 1) * + 4* (F (1 ♦ 1,K- 1) - F(H+1,K+1)) 

GCCA (L, I, J) = .5* (FLOAT ( (1-1) * (J-1) ) * (F (I-J, L- 1) 

1 - F (I- J, L+ 1 ) - F(I*J-2,L-1) + F (I+J-2, L+ 1 ) ) 

2 - FLOAT ( (1-1) * (J+1) ) * (F (I-J-2,L-1) - F (I- J- 2, L+ 1) 

3 - P(I*J,L-1) F { I + J , L ♦ 1 ) ) - FLOAT ( (1 + 1 )* (J- 1) )* (F (I-J+2,L-1) 

4 - F(I-J+2,L + 1) - F(I + J,L-1) + F (1+ J,L+ 1) ) + FLOAT {( 1+ 1 )* (J ♦ 1 ) ) 

5 * ( F (I-J , L- 1 ) - F(I-J,L + 1) - F (I+J + 2,L-1) + F ( 1 + J + 2 , L+ 1 ) ) ) 

GCCB (L , J , I) = -.5* (FLJAT (1-1) **2* (F (I + J-2, L-1) - F (I +J-2 , L+ 1) 

1 + F (I-J, L-1) - F(I-J,L + 1) -F (I + J ,L - 1 ) + F(I+J,L+1) 

2 - F (I- J-2, L- 1 ) + F (I- J-2, L+ 1) ) - FLOAT ( I <■ 1 ) ** 2* 

3 (F (1 + J, L- 1) - F(I + J,L+1) + F (I-J +2 , L- 1 ) - F(I-J+2,L+1) 

4 - F (1 + J + 2, L- 1 ) + ? (1+ J + 2, L+ 1 ) - F (I-J, L-1) + F(I-J,L+1))) 

PI = 3 .14159 

PI 2 = PI*PI 
PI 3 = PI2*PI 
PI 4 = PI3*PI 
RE AD (5 , 1 10 ) AB,NJ 
110 FORMAT (4E10. 3) 

WRITE (6, 1 10 1) A B , .1 U 
1101 FORMAT ( 1 HO , 1P4E10. 3) 

WRITE (6,1) 

READ (5, 120) H V 
WRITE (6, 120) NV 
120 FORM AT ( 15) 

WSITF (6,1) 

1 FORMAT (1H0) 

2 FORMAT ( 1 H 1 ) 

AB 2 = A 3** 2 
AB4 = A 13**4 

C COMPUTE B 

NY = 2*MV+2 
NY = 3 

DO 1C K = 1 , N < 

DO 10 I = 1 , N Y 
L = 2*1 - 1 
DO 1C M = 1 , N X 
DO 10 J = 1 , N Y 


4 
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0040 


N = 2*J - 1 

0041 


AA(K,I,M,J) = I FFF (K,E1) *PF (L,N) + 2.*AB2*FFP (K , M) *FFF(L,N) 



1 + AB4*FF(K,M) *FPFP (L,N) 

0C42 

10 

CONTINUE 


c 

INVFRT AA AND STONE IN AA 

0043 


NAM = NX*N Y 

0044 


1 = 0 

0045 


J = 0 

0046 


DO 512 K = 1, NX 

0047 


DO 512 L = 1 , NY 

0048 


1=1+1 

0049 


DO 508 M = 1 , N X 

0050 


DO 508 N = 1, NY 

0051 


J = J+1 

0052 

508 

AM (I, J) = A A ( K, L, 1, N) 

0053 

512 

J = 0 

0054 


CALL MATIN 2 ( AM , N A 1, 90 , AM W, 90 , A VW , KI NV , 1 8 0) 

0055 


WRITE (6,1) 

0056 


WRITE (6,511) KIN V 

0057 

51 1 

FORMAT (110) 

0058 


1 = 0 

0059 


J = 0 

0060 


DO 522 K = 1 , NX 

0061 


DO 522 L = 1, NY 

0062 


1 = 1 + 1 

0063 


DO 518 M = 1 , NX 

0064 


DO 518 N = 1, NY 

0065 


J = J+1 

0066 

518 

AA(K,L,M,N) = AM (I , J) 

0067 

522 

J = 0 

0068 


DO 14 I = 1 , N V 

0069 


DO 14 K = 1 , N X 

0070 


DO 14 J = 1 , NY 

0071 


L = 2*.t - 1 

0072 


DO 14 M = 1 , N V 

0073 


V(I,K,J,.M) = AB2* (GDC3 (I, K,M) *GCCB (1 , 1 , L) - 2 . *GCC A ( T , K , K) 



1 *GCCA ( 1 ,L, 1 ) + G C U 3 ( I , M,K) * G C C B (1, L, 1) ) *PI4 

0074 

14 

CONTINUE 

0075 


DO 16 K = 1 , N X 

0076 


DO 16 J = 1 , N Y 

0077 


L = 2+ J - 1 

0078 


DO 16 M = 1 , N V 

0079 


DO 16 N = 1 , N V 

0080 


G(K,J,M,N) x 12,* { 1 . - NU**2) *AB2* (GCC A ( K , M , N ) *GCC A (L,1,1) 



1 - GCCB (K , N , M) *GCC J (L, 1 ,1) ) 

0C8 1 

16 

CONTINUE 

0082 


DO 22 I = 1 , N V 

0083 


DO 22 L = I , U V 

0034 


IVSUP1 = I + L 

0085 


DO 22 J = 1 ,N V 

0086 


DO 22 K = 1 , N V 

0087 


IGSUM = J+K 

0088 


B(I,J,K,L) = 0. 

0089 


IS = I+J+K+L 

0090 - 


M = IS-2* (IS/2) 

009 1 


IF(M.NF.O) JO TO 21 
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21 


n ai n 


DATE = 73248 


0092 
009 3 

0094 

0095 

0096 

0097 

0098 

0099 

0100 
0101 
0102 

0103 

0104 

0105 

0106 

0107 

0108 

0109 

0110 
0111 
0112 


IVS = IVSUM-2* (IVSUM/2) + 1 
DO 16 IM = IVS, NX, 2 
IGS = IGSUW-2* (IGS08/2) +1 
DO 18 KK = IGS, NX, 2 
DO 18 JJ = 1 , NY 
DO 18 LL = 1, NY 

18 3(I,J,K,L) = S(I,J,K,L) + V (I , IM , J J , L) * A A (I M , JJ, KK, LL) 

1 *0 (KK, LL, J , K) 

B(I,J,K,L) = -B(I,J,K,L) 

21 CONTINUE 

B(L,J,K,I) = B(I,J,K,L) 

22 CONTINUE 

C WRITE N )NLI NEAR TERMS ONTO TAPE 

WRITE (10) 9 

ENDFILE 10 
WRITF (6,2) 

DO 45 I = 1 , N V 
DO 4 5 J - 1 , N V 
DO 4 5 K = 1 , N V 

45 WRITE (6,47) (B(I,J,K,L), L = 1 , N V ) 

47 FORMAT ( 10E 12. 4) 

STOP 

END 
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DATE 


73248 


000 1 
0002 

0003 

0004 

0005 
000 6 


FUNCTION CC(i<,M) 

CC = 0. 

IF (K . EQ. 3) JC ■= 0J + .5 
IF(K.EQ.-K) c: = JO*. 5 
RETURN . 

END ' 
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fortran 

0001 


0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 
001 1 
0012 

0013 

0014 

0015 

0016 

0017 

0018 
0019 


0020 

0021 

0022 

0023 

0024 

0025 

0026 

0027 

0028 

0029 

0030 

0031 

0032 

003 3 

0034 

0035 

0036 

0037 

0038 

0039 

0040 

004 1 
0042 


IV G LEVEL 21 MATIN2 DATE = 73248 

SUBROUTINE MATIN 2 ( A , N 1 , I A , X, I X, B , INT , N2 ) 

C THIS SUBROUTINE INVERTS THE UPPER LEFT Nl BY N1 CORNER OF 

C MATRIX A , WHICH WHICH HAS AN ACTUAL FIRST DIMENSION OF IA. 

C X AND B ARE DOUBLE PRECISION MATRICES NEEDED FOR WORKING 

C SPACE-X. MUST BE A JO'JBLY DIMENSIONED MATRIX WITH FIRST 

C DIMENSION IX, B IS SINGLY DIMENSIONED AND SHOULD BE OF 

C LENGTH AT LEAST 2 * Nl . INT IS AN INTEGER VARIABLE WHICH 

C IS RETURNED EQUAL TO TWO IF THE MATRIX IS TOO ILL 

C CONDITIONED TO BE INVERTED 

C MODIFIED JORDAN ELIMINATION 

DOUBLE PRECISION 3 (N 2) , X (IX , N 1 ) , PI VOT , T E HP , DA BS 
DIMENSION A (IA, Nl) 

I NT= 1 
N =N 1 

DO 15 1=1, N 
DO 15 J= 1 , N 
15 X (I, J) =A (I, J) 

DO 9 K= 1 , N 

C FIND THE PIVOT 

. PIVOT=0. 

DO 1 I=K , N 
DO 1 J = K , N 

I F(DA3S (X {I, J) ) . Lfi. DABS (PIVOT) ) GO TO 1 
PIVOT=X (I, J) 

A (1,K)=I 
A (2, K) = J 

1 CONTINUE 

'IF(K.EQ.I) CCMP = DA JS (PIVOT) 

IF ( (K. EQ. 1. A ND.COM?. LE. 1 . E-30) .OR. 

1 DABS (PIVOT) . LE. 1. iE-09*COMP) GO TO 14 
C EXCHANGE ROWS 

L = A ( 1 , K ) +1.E-6 
IF(L.EQ.K) GO TO J 
DO 2 J = 1 , N 
TEMP=X (L, J) 

X (L , J ) =X (K , J) 

2 X ( K, J) =TEMP 

C EXCHANGE COLUMNS 

3 L=A (2,K) +1.E-6 
IF(L.EQ.K) GO TO 5 
DO 4 1=1, N 
TEMP=X ( I , L) 

X (I,L)=X (I, K) 

4 X (I , K) =TEMP 

C JORDAN STEP 

5 DO 8 J= 1 ,N 
J2=N+J 

B (J) =1 . D0/PIVOT 
IF (J. NE. K) GO TO 6 
B (J2) =1 . DO 
GO TO 7 

6 B (J) =-X (K, J) *B (J) 

B (J2) =X (J,K) 

7 X (K, J) =0 . 

8 X (J , K) =C . 

DO 9 1=1, N 


8 



FORTRAN IV G LEVEL 21 


HATIN2 


DATE = 73248 


0043 


I2=N+I 


0044 


DO 9 J=1,N 


0045 

9 

X (I, J) “X (I # J) ♦ B (12) *B (J) 


C 

REORDER .FINAL MATRIX 

0046 


DO 1 3 L= 1 , N 


0047 


K = N + 1-L 


0048 


J = A ( 1 , K) *1 . E-6 


0049 


IF(J.EQ.K) GO TO 

1 1 

0050 


DO 10 1=1, N 


0051 


T EM P = X (I, J) 


0052 


X (I , J) =X (I , K) 


0053 

10 

X (I,K) = TfiMP 


0054 

11 

1 = A (2 ,K) +1 . E-6 


0055 


IF(I.EQ.K) GO TO 

1 3 

0056 


DO 12 J»1,W 


0057 


TEP1P=X (I,J) 


0058 


X (I, J)=X(K,J) 


0059 

12 

X (K , J) =TEMP 


0060 

13 

CONTINUE 


0061 


DO 25 1=1, N 


0062 


DO 25 J= 1 , N 


0063 

25 

A (I,J)=X(I,J) 


0064 


RETURN 


0065 

14 

I NT=2 


0066 


RETURN 


0067 


END 
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LEVEL 

21.6 

( MAY 72 > OS/360 FORTRAN H 

COMPILER OPTIONS - NAMB= A AI N , OPT=02 , LI N£CNT=58 r SI2 E=0000K , 

SOURCE, EBCDIC, NOLIST# NO DECK , LOAD, MAP ,NOEDIT,ID, NOIRE P 
C PROGRAM TO COMPUTE AERODYNAMIC ADMITTANCE FUNCTIONS 

C FOR PLATES WITH CLAMPED EDGES 

C EM IS MACH NUMBER 

C DO NOT USE EM = 1.0 

C AB IS PLATE LENGTH/WIDTH RATIO 

C MM AX IS THE NUMBER OF MODES USED IN THE EXPRESSION 

C FOR THE PLATE DEFLECTION 

C IMAX IS THE NUMBER OF POINTS AT WHICH EACH 

C ADMITTANCE FUNCTION IS TO BE COMPUTED 

C THE AEH»S AND AST'S ARE THE ADMITTANCE FUNCTIONS 

ISN 

0002 


DIMENSION ASH (12, 12, 100) ,AEI (12,12,100) 

ISN 

0003 


DIMENSION FEHM (5000) 

ISN 

0004 


READ (5, 1 1 ) EM, AB, M MAX, IMAX 

ISN 

0005 


WRITE (6,11) EM, AB,MNAX, IMAX 

ISN 

0006 

11 

FORMAT (2F10.4, 2110) 

ISN 

0007 


DO 1011 1= 1 , MNAX 

ISN 

0008 


DO 1011 J= 1, FMAX 

ISN 

0009 


DO 1011 K= 1,100 

ISN 

0010 


AEH(I,J,K)- 0. 

ISN 

0011 

1011 

AEI (I,J,K) = 0. 

ISN 

0012 


CIMAX=IM AX 

ISN 

0013 


CMMAX=MMAX 

ISN 

0014 


PI = 3.14159 

ISN 

0015 


SIGF = EM/ ( EM- 1. ) 

ISN 

0016 


X F (EM. GT. 1.) GO TO 15 

ISN 

0018 


EMP = EM* *2/ { 1. -3M**2) 

ISN 

0019 


ABP = (AB**2* 1.) /AB** 2 

ISN 

0020 


SIGF = EiMP * SORT (EMP**2 + EMP*ABP) 

ISN 

0021 

15 

CONTINUE 

ISN 

0022 


DELSIG = SIGF/CIHAX 

ISN 

0023 


WRITE (6,17) SIJF, DELSIG 

ISN 

0024 

17 

FORMAT (2E20. 4) 

ISN 

0025 


GAMMAX=3. *PI 

ISN 

0026 


ALPM AX=SQRT ( (PI* (CMMAX+1.) ) **2+ 100. ) +5. 

ISN 

0027 


DO 24 1=1, IMAX 

ISN 

0028 


CI = I 

ISN 

0029 


S=CI*DSLSIG 

ISN 

0030 


DELG AM=PI/4 . 

ISN 

0031 


DEL ALP=PI/ (4 . * ( 1. +.2*S* (EM+1 . ) /EM) ) 

ISN 

0032 


DEL = -D ELGA M*DELALP/ (PI*EM) **2 

ISN 

0033 


NG AM=G AMM AX/DEL 3 AM 

ISN 

0034 


NALP=ALPMAX/DELALP 

ISN 

0035 


WRITE (6, S 9) N ALP , NGAM 

ISN 

0036 

99 

FORMAT (2120) 

ISN 

0037 


X= DELALP/2. 

ISN 

0038 


DO 22 L= 1 , N A LP 

ISN 

0039 


GERM=0. 

ISN 

0040 


G AM= , 01 

ISN 

0041 


DO 23 K=1 , NGAM 

ISN 

0042 


SQ=SURT (X** 2*SAM**2*A0**2) 

ISN 

0043 


Z=SQ*S/SM 

ISN 

004 4 


CALL GMR (GAM, 1, 1 ,GR,GI) 

ISN 

0045 

45 

C=GR 

ISN 

0046 


GERM=GERM+C*SQ*3J 1 (Z) 
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ISN 

0047 

23 

ISN 

0048 


ISN 

0049 

22 

ISN 

0050 


ISN 

0051 


ISN 

0052 


ISN 

0053 


ISN 

0054 


ISN 

0055 


ISN 

0056 


ISN 

0057 


ISN 

0058 


ISN 

0059 

10 

ISN 

0060 


ISN 

0061 


ISN 

0062 


ISN 

0063 

12 

ISN 

0064 

21 

ISN 

0065 

24 

ISN 

0066 

199 

C 

ISN 

0067 


ISN 

0068 


ISN 

0069 


ISN 

0070 



GAM=GAN't-D2LGA« 

PERM (L) = GJRN 
X=X*DELALP 
DO 21 M=1,MMAX 
DO 21 MH=fl,MNAX 
ALPH= DELALP/2 . 

TERMH=0. 

TEBBI=0. 

DO 10 J-1,NALP 

CALL GMR {ALPH,(1, <R,Gft,GI) 

TERHH=TEHMU<- (GR«SIN (ALPHAS) -GI*COS (ALPHAS)) 

1 *PERM (J) *A1PH 

TERBI=TEBMIf (G8*COS (ALPH*S) ♦GI*SIN (AIPH*S) } *PBRH (J) 

ALPH=ALPH +D2LALP 

AEH (fl,SR,I) -TE8M H*DBL 

AEI (M, MR, I) =TS8MI*DEL 

WHITE (6, 12) AEH(M,Mfl,I) , AEI (M, MR, I) , fl, MR,I 
FORMAT (2 E 20. 3 , 31 1 0) 

CONTINUE 
CONTINUE 
FORMAT (6E20. 3) 

WRITE ADMITTANCE FUNCTIONS ONTO TAPE 
WRITE {10) AEH, AEI, EM, AB, S IGF, IHAX 
ENDFILE 10 
STOP 
END 
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tEVEL 21.6 ( MAY 72 ) 


OS/360 FOHTR AH U 


ISM 

1SN 
ISM 
tS If 
ISM 
ISM 
ISM 
ISM 
ISN 
ISM 
fSM 

ISM 

ISN 

ISM 

ISN 

ISM 

ISM 

ISN 

ISN 

ISN 


COMPILER OPTIONS - N AME= MAIN ,OPT=02, LINECNT=58, SIZE=OOOOK, 


0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 
0011 
0012 

0013 

0015 

0016 

0017 

0018 

0019 

0020 
0021 
0022 


SOURCE, EBCDIC, NOLIST, NODECK, LOAD, MAP, NOEDIT, ID, MOXREF 
SUBROUTINE GMR ( X, H, N , GR,GI) 

C CLAMPED PLATE 

XX = X 

PI = 3. 14159 
AM * M 
AN = N 

A = PI* (A.1-1 , ) 

B = PI*(Ai1+1«) 

C * PI* (AN-1.) 

D = PI* (AN* 1 « ) 

14 CONTINUE 

DEKOM = (X**2-A**2) * (X**2-B**2) *(X**2-C**2) 

1 * (K**2-l)**2) 

IF (AB3 (OENOM) .LT. 1 .0S-10) GOTO 12 

GR = AMP*( (1 . * (- 1.) ** (M*N) ) ♦ ( (-1. ) **H* (- 1. )**M) 

1 *COS(X)) 

GI = A MP*5I N (X) * ( (-1.) **N - (-1.)**M) 

X = XX 
RETURN 

12 CONTINUE 

X = X+.01 
GO TO 14 
END 


12 



tRVEl 21.6 ( HAY 72 ) 


OS/360 fORTBAN U 


ISN 0002 

ISH 0003 
ISN 0004 
ISH 0005 

ISN 0006 
ISN 0007 
ISN 0008 

ISN 0009 


COMPILER OPTIONS - NAME= MAI N ,OPT*02 , LINECNT*58, SIZE*OOOOK, 

30 JRCE, EBCDIC, NOLIST, NODECK, LOAD, HAP # NOEDIT, ID, MOXBIF 
FUNCTION BJ 1 (X) 

C POLYNOMIAL APPBOXIHATION POB BESSEL PONCTION 

IF (X-3.) 1,2,2 

1 Y=(X/3.)**2 

BJ 1 =X* (.5-. 56249 985* I*. 2 1093573*T**2«. 03954289*T**3* 

1 . 00443 31 9*T* *4-. 0 003 17 61 *Y**5+. 00001 109*Y*<6) 

GO TO 3 

2 T-3./X 

F1-. 797884S6+.0GJ001 56*Y*. 01 65966 7*Y**2*. CO0 17 105* Y** 3- 
1 .00249S1 1*Y**4*. JO 1 1 36 5 3* Y**5- . 00020033*T**6 

THl=X-2. 356 1 9 44 4+ . 12499612+Y+. 0C0 0S656*Y**2-. 006 3787 9*Y**3* 


ISN 0010 
ISN 0011 
ISN 0012 
ISN 0013 


. 0 0074 348* Y**4+. 0 007 9824* Y**5-. 00 029166*1**6 
BJ1=F1*COS (TH1) /3QBT (X) 

CONTINUE 

RETURN 

END 
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LEVEL 

21.6 

( MAY 72 

> 

OS/36Q FORTRAN H 


COMPILER 

OPTIONS - NAM & = MA I N ,0 PT=02 , LI NECNT=58 ,SI ZE=OOOOK , 





SOURCE , EBCDIC, NOLIST, NODECK, LOAD, SAP, NOEDIT , ID, NOXREP 



c 


FLUTTER PROGRAM FOR CLAMPED-EDGE PLATES USING 



c 


LINEARIZED POTENTIAL PLOW AERODYNAMICS 



c 


LAMDA=DYN AMIC PHSSSUB 3, MU=FLOW DENSITY, MACH=NACH NUMBER 



c 


AB=PLATE LENGTH/WIDTH RATIO, RXA, RYA=APPLI£D IN-PLANE 



c 


LOAD (POSITIVE IN TENSION) , PSTAT=STATIC PRESSURE 



c 


DIFFERENTIAL, CAVITY=CAVITY ACOUSTIC PARAMETER 



c 


DAM P= STRUCTURAL DAMPING FACTOR 



c 


NV = #OF MODES, H = I INTEGRATION STEP INTERVAL, TPBINT=PRINT-OUT 



c 


INTERVAL, TFINAL=TIME AT WHICH INTEGRATION STOPS 



c 


SCALE= MAXIMUM ANTICIPATED DEFLECTION <POR GRAPH ROUTINE) 



c 


X, , . = ALPHANUMERIC CHARACTERS FOR GRAPH ROUTINE 



c 


THE A * S ARE MODAL AMPLITUDES 



c 


THE W'S ARE THE PANEL DEFLECTION AT 15 EVENLY SPACED POINTS 



c 

c 


ALONG THE PANEL CENTERLINE 



c 

c 


REMOVE CAaDS #177 THROUGH 189 FOR ZERO EDGE RESTRAINT 



c 

c 


CALCULATION 

ISN 

0002 

c 


REAL LAMDA , MU 

ISN 

0003 



REAL NU 

ISN 

0004 



REAL MACH 

ISN 

0005 



DIMENSION B ( 12, 12, 12, 12) 

ISN 

0006 



DIMENSION ASH (12, 12, 10 0) , A El (12, 1 2, 100) 

ISN 

0007 



DIMENSION AS (500, 12) , DAS (500, 12) 

ISN 

0008 



DIMENSION S ( 12, 12) ,C (12, 12) ,D (12, 12) ,PHIX(12, 12) ,PHIY(12, 12) 

ISN 

0009 



DIMENSION A ( 12) , DA (12) , DDA ( 1 2) , DD AS (4 , 1 2) 

ISN 

0010 



DIMENSION 3(12, 12) 

ISN 

0011 



DIMENSION W (15) , F (12) 

ISN 

0012 



DIMENSION 3T (12, 12) 

ISN 

0013 



DIMENSION MM (12, 12) 

ISN 

0014 



DIMENSION W V (24) 

ISN 

0015 



DIMENSION RI NE (6 1 ) 

ISN 

0016 



REAL*8 WM 

ISN 

0017 



REAL*8 ktf 

ISN 

0018 



PP (I , M) = CC (I- 1 , M-1) -CC (1-1 ,M + 1) -CC(I+1,M-1) *CC{I+1,M*1) 

ISN 

0019 



PPX(I,J) * -PI* FLOAT (J- 1) * (CS (I- 1 , J - 1 ) -CS (1 + 1, J-1) ) 




1 

♦ PI*FLJAT (J + 1) * (CS (1-1 ,J+ 1) -CS (I+1,J* 1) ) 

ISN 

0020 



PPXX (I , M) - -?I2*FLOAT (B-1) **2* (CC (I- 1,M-1) -CC (1*1, M-1) ) 




1 

♦ PI2*F„OAT { M+ 1 ) **2 * (CC (1-1 , M + 1 ) -CC (I* 1,M* 1) ) 

ISN 

0021 



PPXXXX (I , M) - PI4* FLOAT (M-1) ** 4* (CC (I - 1, M-1) -CC (1 + 1 ,H-1) ) 




1 

- PI4*FL0AT (M*1) **4* (CC (I- 1 , M+ 1 ) - CC (1+ 1 , M* 1 ) ) 

ISN 

0022 

1 


FORMAT ( 1 H 0 ) 

ISN 

0023 

2 


FORMAT ( 1 H 1 ) 

ISN 

0024 



PI = 3. 14 159 

ISN 

0025 



PI2 = PI*P1 

ISN 

0026 



PI 3 = PI 2* PI 

ISN 

0027 



PI 4 = PIJ*PI 

ISN 

0028 



NU = . 3 

ISN 

0029 



READ (5,701) CROSS, BLANK, DOT, SCALE 

ISN 

0030 

701 


FORMAT ( 3 A 1 , F7.1) 

ISN 

0031 



WRITE (o, 700 ) CROSS , BLANK, DOT, SCALE 

ISN 

0032 

700 


FORMAT (K, 3A1, F7.2) 
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ISN 0033 


SCALE = 30. /SCALE 

ISN 0034 


DO 702 I = 1,61 

ISN 0035 

702 

BINE (I) - ELANK 


C 

READ NONLINEAR TERMS FROM TAPE 

ISN 0036 


READ (10) 3 

ISN 0037 


REWIND 10 


c 

READ AERODYNAMIC ADMITTANCE FUNCTIONS 

ISN 0038 


READ (12) AEH,AEI, MACU, AS, SIGF, ISM AX 

ISN 0039 


REWIND 12 

ISN 0040 


WRITE (6,1) 

ISN 0041 


WRITE (6,13) MACS, AB,N V, IS BAX, SIGF 

ISN 0042 

13 

FORM AT (2F1 0.4, 21 10, El 0.3) 

ISN 0043 

5 

FORMAT (1P8E9.2) 

ISN 0044 

8 

FORMAT (10E1.2.3) 

ISN 0045 


WRITE (6,1) 

ISN 0046 


READ (5,110) LAMDA, MU, PSTAT, CAVITY 

ISN 0047 

110 

FORMAT (4E10. 3) 

ISN 0048 


WRITE (6,1101) LAMDA, MU, PSTAT, CAVITY 

ISN 0049 

1101 

FORMAT ( IdO, 1P4E10. 3) 

ISN 0050 


READ (5,110) BX A, RY A 

ISN 0051 


WRITE (6, 1 101) R X A , R Y A 

ISN 0052 


WRITE (6,1) 

ISN 0053 


READ (5,116) DAMP 

ISN 0054 

116 

PORHAT (£13. 3) 

ISN 0055 


WRITE (6,117) DAMP 

ISN 0056 

117 

FORM AT ( ' STRUCTURAL DANPING= ',E12.3) 

ISK 0057 


WRITE (6,1) 

ISN 0058 


READ (5,1113) MV, H,TPRINT, TFINAL 

ISN 0059 


WRITE (6,1113) NV,H,TP2INT,TPINAL 

ISN 0060 

1113 

FORMAT (I10,3E1J. 3) 

tSN 0061 


WRITE (6,1) 

tSN 0062 


READ (5,113) ( A ( I ) , I = 1 , NV ) 

ISN 0063 


WRITE (6,115) (A (I) , I = 1 , N V) 

ISN 0064 


READ (5,113) (DA (I) , I = 1,NV) 

TSN 0065 


WRITE (6,115) ( 3 A ( £) , 1= 1,NV) 

ISN 0066 

113 

FORMAT (6 E 1 0 , 3) 

ISN 0067 

115 

FORMAT ( 1X,6E12. 3) 

ISN 0068 


AB2 = A 3**2 

ISN 0069 


AB4 = Ad**4 

ISN 0070 


RISMAX = ISMAX 

ISN 0071 


DELS = 5IGF/RISMAX 

ISN 0072 


ROOT = SORT (LAMJA/MU) 

ISN 0073 


HAERO = RGOT*H 

ISN 0074 


NAERO = DELS/HA 3RO 

ISN 0075 


IF (NAERO. ul. 1) GO TO 404 

ISN 0077 


IF (NAERO. GT. 20) NAERO = 20 

ISN 0079 


DELSIG = N AESO*<3 AERO 

ISN 0080 


IM AX = S IGF/DELS IG 

ISN 0081 


GO TO 406 

ISN 0082 

404 

CONTINUE 

ISN 0083 


NAERO = 1 

ISN 0084 


DELSIG = DELS 

ISN 0085 


H = DELSIG/RCOT 

ISN 0086 


IMAX = ISMAX 

ISN 0087 

406 

CONTINUE 

ISN 0083 


HP = D2LSIG/ROOT 


PROM TAPE 


15 - 



ISN 0089 


WHITE (6,1) 

ISN 0090 


WRITE (6,4 31) H,DELSIG,NAERO,IMAX 

ISN 0091 

401 

FORMAT (2H20. 3,215) 

ISN 0092 


IF {IMAX.GT. 100) STOP 

ISN 0094 


NSTORE = NABflO*IMAX 

ISN 0095 


DO 3 1= 1 , N V 

ISN 0096 


DO 3 J= 1 , N V 

ISN 0097 


DO 3 K- 1 , ISN AX 

ISN 0098 


AEH ( I, J, 100-K+1) = AEH (I, J,ISMAX-K+1) 

ISN 0099 


AEI (I,J, 1C 0-K ♦ 1 ) = AEI (I,J,ISHAX-K+1) 

ISN 0100 

3 

CONTINUE 

ISN 0101 


IP - 1 00-ISM AX+ 1 

ISN 0102 


DO 400 M = 1 , N V 

ISN 0103 


DO 400 N = il,NV 

ISN 0104 


DO 400 I = 1 , INAX 

ISN 0105 


X = PLOAT (I) *DELSIG/DELS 

ISN 0106 


J = INT(X) 

ISN 0107 


P=X-AINT ( X) 

ISN 0108 


JP=IP-1+J 

ISN 0109 


IF (J) JOO,3CO,JQ1 

ISN 0110 

300 

AEH (M,M,I) =AEH(1, N,JP*1) *P 

ISN 0111 


AEI (M,N,I)=AEI(.M,N,JP*1) * P 

ISN 0112 


GO TO 400 

ISN 0113 

30 1 

AEH (N, N ,1) = A3d(M,N,JP)*{1.-P) + AEH ( N, N , JP+ 1 ) *P 

ISN 0114 


AEI(M,N,I> = A8L(H, N,JP)* (1.-P) «■ AEI (M, N, JP+1) *P 

ISN 0115 

400 

CONTINUE 

ISN 0116 


DO 410 A = 1 , NV 

ISN 0117 


DO 410 A = M,NV 

ISN 0118 


DO 4 10 I = 1 , III AX 

ISN 0119 


AEH(N,M,I) - (-1 . ) ** {M+N) *AEH (N,N,I) 

ISN 0120 


AEI (N,n,I) = (- 1. )** (H + N) *AEI(H,M,I) 

ISN 0121 

410 

CONTINUE 

ISN 0122 


DO 30 I = 1,8V 

ISN 0123 


DO 30 J -= 1 , NV 

ISN 0124 


S(I,J) - PP(I,J) *PP(1,1) 

ISN 0125 


ST (I, J) = S(I,J) 

ISN 0126 


C(I,J) - PPXXXXfl, J) *PP (1,1) ♦ 2.*AB2*PPXX(I,J)*PPXX(1 


1 

♦ AB4*PP (I, J) *PPXXXX (1, 1) 

ISN 0127 


D (I, J) = PPX (I, J) *PP (1,1) 

ISN 0128 


PHI X (I , J) =-PPXX (I,J) *PP(1,1) 

ISN 0129 


PHI I (I,J) =-PP(I, J) *PPXX (1,1) 

ISN 0130 

30 

CONTINUE 

ISN 0131 


WRITE (6,1) 

ISN 0132 


DO 910 I - 1 , NV 

ISN 0133 

910 

WRITE (6,5 7 3) ( S(I,J), J = 1,NV) 

ISN 0134 


WRITE (6,1) 

ISN 0135 


DO 580 I = 1 , N V 

ISN 0136 

530 

WRITE (6,57 3) (C(I,J), J = 1,NV) 

ISN 0137 

573 

FORMAT (8E16.4) 


C 

INVEST 3 AND STORE IN S 

ISN 0138 


CALL MATIN2 (S,NV, 12, N M , 1 2 , W V , KI NV RT, 24) 

ISN 0139 


WRITE (6, 1) 

ISN 0140 


WRITS (6, 5067) KINVRT 

ISN 0141 

5067 

FORMAT (15) 

ISN 0142 


WRITE (6,1) 

ISN 0143 


DO 912 I = 1 , N V 


16 



ISN 0144 

912 

WRITE(6,57i) ( D(I,J), J * 1,NV) 

ISN 0145 


WRITE (6, 1 ) 

ISN 0146 


DO 914 I = 1 , N V 

ISN 0147 

914 

WRITE (6,57 3) (PHIX(I,J>, J = 1,NV) 

ISN 0148 


WRITE <6, 1) 

ISN 0149 


DO 916 I = 1,N7 

ISN 0150 

916 

WRITE (6, 573) (P8IY(I,J), J= 1,NV) 

ISN 0151 


WRITE (6,1) 


C 

SET UP INTEGRATION 

ISN 0152 


T = 0. 

ISN 015.1 


TP •= 0. 

ISN 0154 


DO 70 I - 1,500 

ISN 0155 


DO 70 J = 1, NV 

ISN 0156 


AS (I , J) = 0, 

ISN 0157 

70 

DAS (I,J) = 0. 

ISN 0158 


DO 74 I * 1,4 

ISN 0159 


DO 74 J = 1, NV 

ISN 0160 

74 

DD AS (I,J) = 0. 

ISN 0161 


IMAXP = INAX - 1 

ISN 0162 


HHH = H/24. 

ISN 0163 


WRITE (6,2) 


C 

PREDICTOR ROUTINE 

ISN 0164 

56 

CONTINUE 

ISN 0165 


IF (A (1) ,GT. 10.) STOP 


C 

COMPUTE SECOND DIRIVATIVES 

ISN 0167 


DO 490 I = 1 , N V 

ISN 0168 


DO 490 J - 1 ,NV 

ISN 0169 


Q (1,0) * (D(I,J)*A(J) 4- ST (I, J) *DA (J) /ROOT) /MACH 

ISN 0170 


DO 488 K - 1 , 1 .1 A X P 

ISN 0171 


KP = NA£R3*K 

ISN 0172 


Q (1,0) — Q (I , J) +-DA3 (KP,J) *AEI (J, I , K) *HP | 



1 4-AS (KP,J) *A£H (J, I,K) *DELSIG I 

ISN 0173 

488 

CONTINUE 

ISN 0174 

490 

CONTINUE 

ISN 0175 


SX = 0. 

ISN 0176 


SY = 0. 


C 

REMOVE THR FOLLOWING 12 CARDS FOR A ZERO EDGE RESTRAINT 


C 

CALCULATION 

ISN 0177 


SY = . 5 + A ( 1) **2 

ISN 0178 


DO 212 fl = 1 , N V 

ISN 0179 


rh = n 

ISN 0180 


SX = SX ♦ (RM**2<-1.) *A (H) **2 

ISN 0181 

212 

SY - SY ♦ A (il) **2 

ISN 0182 


N VP = NV-2 

ISN 0183 


IF (NVP . LT. 1) GO TO 21 5 

ISN 0185 


DO 214 N =■ 1 , N VP 

ISN 0186 


RM = PI 

ISN 0187 


SX - SX - (RM+1.) **2*A(N)*A(fl*2) 

ISN 0188 

214 

SY = S Y - A ( 8) * A (114-2) 

ISN 0189 

215 

CONTINUE 

ISN 0190 


AB2RXB = 12. *PI2* (.75*SX ♦ NU*AB2*ST) 4- RXA 

ISN 0191 


RYB = 12. *PI2* (A32*SY ♦ NU*.75*SX) ♦ RYA 

ISN 0192 


DO 250 I = 1 , NV 

ISN 0193 


F(I) = 0. 

ISN 0194 


DO 216 M = 1 , N V 

ISN 0195 

216 

F(I) = F (I) - AB2RXB*PHIX(I,M) * A ( H) - AB2*RYB*PHIY (I , (1) *A (M) 
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ISN 0196 


DO 220 J = 1,NV 

ISM 0197 


DO 220 K = 1,NV 

ISN 0198 


LSUN = I + JtK 

ISN 0199 


LS = 2* (LSUil/2) -LSU1U2 

ISN 0200 


DO 220 L = LS,NV,2 

ISN 0201 

220 

F(I) = F (I) - B(I,J,K r L)*A(J)*A(K)*A(L) 

ISN 0202 


DO 230 J - 1,N 7 

ISN 0203 

230 

F (I) = F(I) - C(I,J)*A(J) - LAMDA*Q(I,J)-DAHP*C(I,J) *DA(J) 

ISN 0204 

250 

CONTINUE 

ISN 0205 


F ( 1 ) = F ( 1 ) - CAVITY* A ( 1) ♦ PSTAT 

ISN 0206 


DO 240 I - 1 » N V 

ISN 0207 


DDA(I) = 0. 

ISN 0208 


DO 240 J = 1 , N V 

ISN 0209 


DDA(I) = DDA(I) + S(L,J)*F(J) 

ISN 0210 

240 

CONTINUE 


C 

PRINT OUTPUT 

ISN 0211 


IF (T, LT. TP) GO TO 350 

ISN 0213 


TP = TP + TP SI NT 


c 

PRINT NODAL AMPLITUDES, VELOCITIES, AND ACCELERATIONS 

ISN 0214 


WRITE (6,345) T 

ISN 0215 

345 

FORMAT (6H IIM2=,F7.4) 

ISN 0216 


WRITE (6, 347) (A (I), I - 1*HV) 

ISN 0217 


WRITE (6 , 347) (JA (I) , I * 1,NV) 

ISN 0218 


WRITE (6 , 347) (LJA(I), I = 1,NV) 

ISN 0219 

347 

FORMAT (6E1 1. 3) 

ISN 0220 


DO 348 I = 1,15 

ISN 0221 


W(I) = 0. 

ISN 0222 


THETA = PI *F LOAT { I) / 1 5 • 

ISN 0223 


DO 346 J = 1 ,NV 

ISN 0224 

346 

W(I) * W(I) ♦ 2. ♦A (J) * (COS (FLOAT (J-1) *THETA) - COS ( FLOAT (J *■ 1 ) 



1 *THETA)> 

ISN 0225 

348 

CONTINUE 


C 

PRINT PLATE DEFLECTION AT 15 EQUALLY SPACED POINTS ALONG THE 


C 

CENTERLINE OF THE PANEL 

ISN 0226 


WRITE (6,349) (i ( I) r 1-1 ,15) 

ISN 0227 

349 

FORNAT (8F7.2) 

ISN 0228 


R IN E ( 1 ) = DOT 

ISN 0229 


RINE (31) = DOT 

ISN 0230 


RINE(fil) ^ DOT 

ISN 02 31 


L = SCALE* W (12) 

ISN 0232 


LP - 3HL 

ISN 0233 


IF (IASS (L) » L E. 33) RINE (LP) = CROSS 


C 

GRAPH DEFLECTION OF POINT ON LATERAL CENTERLINE OF PANEL 


C 

3/4 OF WAY EACH FROM LEADING EDGE 

ISN 0235 


WHITE (6,703) (EINB(I), I = 1,61) 

ISN 0236 

703 

FORNAT (68X, 61A1) 

ISN 0237 


IF (IABS (L) .LE. 33) HINS(LP) = BLANK 

ISN 0239 


IF (T.GE, TFiNAL) GO TO 57 

ISN 0241 

350 

CONTINUE 


C 

STORE VARIABLES 

ISN 0242 


DO 24 J = 2, NST ORE 

ISN 0243 


K = NSTORE-J *2 

ISN 0244 


KP = K-1 

ISN 0245 


DO 24 I = 1,NV 

ISN 0246 


AS (K,I) * AS (KP, I) 

ISN 0247 


DAS (K, I) = DAS (KP, I) 
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XSN 0248 
ISN 0249 
ISN 0250 
ISN 0251 
ISN 0252 
ISN 0253 
ISN 0254 
ISN 0255 
ISN 0256 
ISN 0257 
ISN 0258 
ISN 0259 

ISN 0260 
ISN 0261 

ISN 0262 

ISN 0263 
ISN 0264 
ISN 0265 
ISN 0266 
ISN 0267 
ISN 0268 


24 CONTINUE 

DO 26 J = 2,4 

K = 6-J 

KP * K- 1 

DO 26 I = 1 , N V 

DDAS (K, I) = DDAS ( KP, I ) 

26 CONTI N US 

DO 28 I - 1 , N V 
AS (1,1) = A ( I) 

DAS (1,1) = D A (I) 

DDAS (1 , I) * IDA ( I) 

28 CONTINUE 

C PHEDICT 

DO 20 I - 1, NV 

A (I) * A(I) ♦! HH* (55.*DAS (1 ,1) -59. *DAS (2,1) 

1 *37.*DAS(J,I)-9.*DAS (4,1)) 

DA (I ) = DA (I) fHUH* (55.*DDAS ( 1 , 1) - 59« *DDAS <2 # I) 
1 *37.*DDAS (3, I) -9. *DOAS (4,1)) 

20 CONTINUE 

T * T+H 
GO TO 56 
57 CONTINUE 

STOP 
END 


<3 
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OS/360 FORTRAN H 


,£VEL 

21.6 ( 

MAY 72 

) OS/360 FORTRAN H 


COMPILER 

OPTIONS - NASE= MAlN,OPr=02,LINECNT=58,SIZB=0000N, 




SO a RC2, EBCDIC, NOLIST , NODECK, LOAD, BAP, NOEDIT, ID, NOXREF 

ISN 

0002 


SUBROUTINE MATIN 2 (A, Ml , IA, X, IX, B # INT, N2) 



c 

THIS SUBROUTINE INVERTS THE UPPER LEFT N1 BY N1 CORNER OF 



c 

MATRIX A, rfHICH 2HICH HAS AN ACTUAL FIRST DIMENSION OF IA. 



c 

X AND B ARE DOUBLE PRECISION MATRICES NEEDED FOR BORNING 



c 

SPACB-X MUST BE A DOUBLY DIMENSIONED MATRIX ilTH FIRST 



c 

DIMENSION IX, B IS SINGLY DIMENSIONED AND SHOULD BE OF 



c 

LENGTH AT LEAST 2*K1. INT IS AN INTEGER VARIABLE NHICH 



c 

IS RETURNED EQUAL TO T MO IP THE MATRIX IS TOO ILL 



c 

CONDITIONED TO BE INVERTED 



c 

MODIFIED JORDAN ELIMINATION 

ISN 

0003 


DOUBLE PRECISION B (N2) , X (IX , N 1 ) , PI VOT, TEMP, DABS 

ISN 

0004 


DIMENSION A (I A, Nl) 

ISN 

0005 


INT* 1 

ISN 

0006 


N = N1 

ISN 

0007 


DO 15 1=1, N 

ISN 

0008 


DO 15 J=1 , N 

ISN 

0009 

15 

X(I,J)=A(I,J) 

ISN 

0010 


DO 9 K= 1 , N 



c 

PIND THE PIVOT 

ISN 

0011 


PI VOT=0. 

ISN 

0012 


DO 1 I=K, N 

ISN 

0013 


DO 1 J=K , N 

ISN 

0014 


IF (DABS (X (I, J) ) .LB. DABS (PIVOT)) GO TO 1 

ISN 

0016 


PIVOT=X (I,J) 

ISN 

0017 


A(1,K)=I 

ISN 

0018 


A(2,K)=J 

ISN 

0019 

1 

CONTINUE 

ISN 

0020 


IF(K.EQ.I) COMP-DABS (PIVOT) 

ISN 

0022 


IF( (K.EQ. 1. AND.COMP.LB. 1.E-30) .OR. 




1 DABS (PIVOT) . LE. 1 . 5E-Q9+COMP) GOTO 14 



c . 

EXCHANGE HONS 

ISN 

0024 


L=A ( 1 , K ) fl.E-6 

ISN 

0025 


IF(L.EQ.K) GO TO J 

ISN 

0027 


DO 2 J= 1 , N 

ISN 

0028 


TEMP=X (L, JJ 

ISN 

0029 


X (L, J) =X (K, J) 

ISN 

0030 

2 

X (K, J) =TEMP 



C 

EXCHANGE COLUMNS 

ISN 

0031 

■3 

L= A (2 , K) ♦ 1 . E-6 

ISN 

0032 


IF(L.EQ.K) GO TO 5 

ISN 

0034 


DO 4 1=1, N 

ISN 

0035 


TEMP=X (I,L) 

ISN 

00.16 


X(I,L)=X(I,K) 

ISN 

0037 

4 

X (I,K) =TSBP 



C 

JORDAN STEP 

ISN 

0038 

5 

DO 8 J=1,N 

ISN 

00 39 


J2=N*J 

ISN 

0040 


B(J)=1. D0/PIVOT 

ISN 

0041 


IF(J.NE.K) GO TO 6 

ISN 

0043 


B ( J 2) = 1 . DO 

ISN 

0044 


GO TO 7 

ISN 

0045 

6 

B(J) =-X (K,J) * E (J) 

ISN 

0046 


B ( J 2 ) =X (J,K) 

ISN 

0047 

7 

X (K , J) =0. 
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ISN 

0048 

8 

X (J,K)=0. 


ISN 

0049 


DO 9 1=1, N 


ISN 

0050 


I2=N+I 


ISN 

0051 


DO 9 J=1,N 


ISN 

0052 

9 

X(I,J)=X(I,J> ♦B{I2) *B<J) 



C 

REORDER PI HAL HATHIX 

ISN 

0053 


DO 13 L= 1 , N 


ISN 

0054 


K=N ♦ 1 - L 


ISN 

0055 


J= A (1,K) M.E-6 


ISN 

0056 


IF(J.EQ.K) GO TO 

1 1 

ISN 

0058 


DO 10 1=1, N 


ISN 

0059 


TEH P=X (I, J) 


ISN 

0060 


X (t,J)=X (I,K) 


ISN 

0061 

10 

X (I,K) = T8HP 


ISN 

0062 

11 

1 = A (2, K) +1.E-6 


ISN 

0063 


IP (I. EQ. K) GO TO 

13 

ISN 

0065 


DO 12 J=1,H 


ISN 

0066 


TEHP=X (I,J) 


ISN 

0067 


X (I , J) = X (K, J) 


ISN 

0068 

12 

X (K,J)=TEHP 


ISN 

0069 

13 

COHTINOE 


ISN 

0070 


DO 25 1=1, H 


ISN 

0071 


DO 25 J= 1 , N 


ISN 

0072 

25 

A (I, J) =X (I, J) 


ISN 

0073 


RETURN 


ISN 

0074 

14 

INT=2 


ISN 

0075 


RETURN 


ISN 

0076 


END 
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OS/360 FORTRAN 8 


LEVEL 21.6 ( MAY 72 ) 

COMPILER OPTIONS - NAME- MAIN, OPT=Q2, LINECIT*58 ,SIZI*OQQOK, 

SOURCE, EBCDIC, NOLIST, NODICK, LOAD, HAP, NOEDIT, ID, NOXBEF 
ISN 0002 FONCTION CS(M,N) 

ISN 0003 CS = 0. 

ISN 0004 PI = 3. 14159 

ISN 0005 RM = M 

ISN 0006 RN = N 

ISN 0007 IF(MOD{M+N,2) . N iJ. 0) CS * -2 . *IN/ (PI* (RH**2“RN**2) ) 

ISN 0009 RETURN 

ISN 0010 END 


22 



LEVEL 21.6 ( BAT 72 ) 


OS/360 FORTRAN H 


ISN 0002 
ISM 0003 
ISN 0004 
ISN 0006 
ISN 0003 
ISN 0009 


COMPILER OPTIONS - NANE* BAH ,OPT=02, 1INECNT*53,SIZE*0000K, 

SOURCE, EBCDIC, ROUST, ROD SCK, LOAD , MAP, NOEDIT, ID, NOXREF 
FUNCTION CC (K , N) 

CC = 0. 

IF(K.EQ.fl) CC = CO. 5 
IF(K.EQ.-fl) CC * CO. 5 
RETURN 
END 



LEVEL 

21.6 

( BAT 72 

) OS/360 FORTRAN H 


COMPILER 

OPTIONS - NAME* MA IN , OPT=02 , LINECNT*58, SIZB*OOOOK, 




SOURCE, EBCDIC, NOLIST, NODECK, LOAD, HAP, NOBDIT, ID, NOXHBF 



c 

PROGRAM TO CCHPUTB STRESSES DOB TO PLATB DEFLECTION FOR 



c 

CLAMPED PLATS MI TH COMPLETE IN-PLANE EDGE RESTRAINT 



c 

A B* PLATE LENGTH/WIDTH RATIO, NY*# OF BODES, NU*POISSON * S RATIO 



c 

XA, YB*COORi)INATES OF POINT AT STRESSES ARE CONFUTED 



c 

(XA, IB ARE HONDIMENSIOKALIZBD BT PANEL LENGTH 6 WIDTH) 



c 

XA=H/1 6,N=Q, 1,2,**. 17 



c 

TB IS SPECIFIED AS INPUT DATA 



c 

THE A*S ARE THE MODAL AMPLITUDES (WHICH DEFINE THE PLATE 



c 

DEFLECTION) 

ISN 

0002 


REAL NU 

ISN 

0003 


DIMENSION A (12) 

I SN 

0004 


DENOHI(K) = (FLOAT (K* *2) ♦ 16.*AB2)**2 

ISN 

0005 


DENOM2 (K) * (FLOAT (K**2) ♦ 4,*AB2)**2 

ISN 

0006 


B1 (N,N) = -FLOAT (2*M* (M+-N) ♦ 4) /DENOM1 (M-N) 

ISN 

0007 


B2(M,N) = FLOAT ( (M-1) * (M+N) )/DBNOMl (M-N-2) 

ISN 

0008 


B3(M,N) = FLOAT ( (N+ 1 ) * (Mt-N) ) /DENOH1 (M-N* 2) 

ISN 

0009 


B4(H,N) * FLOAT (2*M* (M-N) ♦ 4) /DENON1 (fl*H) 

ISN 

0010 


B5(M,N) = -FLOAT ( (M-1) * (M-N)) /DENOH1 (H4-N-2) 

ISN 

0011 


B6(H,K) = -FLOAT ( (M + 1 ) * (M-N)) /DENOM1 (H+N+2) 

ISN 

0012 


B7(H,N) = FLOAT (4*H**2 * 4) /DENOH2 (M-N) 

ISN 

0013 


B8(M,N) = -FLOAT (2* (M-1) **2)/DENOM2 (H-N-2) 

ISN 

0014 


B9(H,N) * -FLOAT (2* (M*1) ** 2 ) /DENOM2 (M-N >2) 

ISN 

0015 


BIO (M,N) = -FLOAT (4*M**2 ♦ 4) /DENOM2 (M+N) 

ISN 

0016 


B 1 1 (M,N) * FLOAT (2* ( M- 1 ) **2) /DENOH2 (M* N- 2) 

ISN 

0017 


B1 2 (H, N) = FLOAT (2*(M*1)**2)/DENOM2(M*N*2) 

ISN 

0018 


CS (M, X) = COS (FLOAT (M) *PI*X) 

ISN 

0019 


READ (5,697) AB, NV 

ISN 

0020 


WRITE (6,697) A3, NY 

ISN 

0021 

697 

FORMAT (P 10. 3, 110) 

ISN 

0022 


READ (5,697) IB 

ISN 

0023 


WRITE (6, 6S7) IB 

ISN 

0024 


READ (5,699) (A (I) , I * 1,NV) 

ISN 

0025 


WRITE (6,699) (A (I), I = 1,NV) 

ISN 

0026 

699 

PORHAT (6F10.4) 

ISN 

0027 


PI * 3.14159 

ISN 

0028 


PI2 = PI**2 

ISN 

0029 


PI4 = PI**4 

ISN 

0030 


AB2 = AB**2 

ISN 

0031 


AB4 = AB**4 

ISN 

0032 


NU = .3 

ISN 

0033 


XA = 0. 

ISN 

0034 


CHI = 1. - CS (2,TB) 

ISN 

0035 


DO 710 II = 1,17 

ISN 

00 36 


WRITE (6,1) 

ISN 

0037 

1 

FORMAT (ISO) 

ISN 

0038 


WRITE (6,701) XA, IB 

ISN 

0039 

701 

FORMAT (2320. 3) 

ISN 

0040 


W * 0. 

ISN 

0041 


WXX = 0. 

ISN 

0042 


WYY = 0. 

ISN 

0043 


DO 704 M = 1 , NV 

ISN 

0044 


PSI * CS (il- 1 , XA) - CS (M+ 1 , X A) 

ISN 

0045 


RH = M 

ISN 

0046 


PSIXX * - (RN-1.) **2*PI2*CS (H-1,XA) + (R M* 1 . ) **2 *PI2* 
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704 


ISN 0047 
ISN 0048 
ISM 0049 
ISN 0050 
ISN 0051 
ISN 0052 
ISN 0053 
ISN 0054 
ISN 0055 
ISN 0056 

ISN 0057 
ISN 0058 
ISN 0059 
ISN 0060 
ISN 0061 
ISN 0062 
ISN 0063 
ISN 0064 
ISN 0066 
ISN 0067 
ISN 0068 
ISN 0069 
ISN 0070 
ISN 0071 
ISN 0072 
ISN 0073 
ISN 0074 
ISN 0075 
ISN 0076 
ISN 0077 
ISN 0078 


ISN 0079 
ISN 0080 


ISN 0081 
ISN 0082 
ISN 0083 
ISN 0084 


ISN 0085 
ISN 0086 


ISN 0087 
ISN 0088 


706 


708 

712 


1 CS(M+1,XA) 

W = W + A (M) *PSI 

WXX = WXX + PSIXX* A (tl) 

W YY = WYY + PSI*A(M) 

W = W*CHI 

WXX = WXX*CHI 

WYY = 4.*PI2*CS(2,YB) *WYY 

WRITE (6,701) W 

WRITE (6,701) WXX, WYY 

SIGHAX = -, 5* (WXX + N SJ* AB2*WY Y) 

SIGH AY = -.5* (A32*WY Y *■ NU*WXX) 

PRINT PLATE BENDING STRESSES AT POINT (XA,YB) 

WRITE (6,701) SIGHAX, SIGMAY 
SX = 0, 

SY = . 5*A ( 1 ) * *2 
DO 706 M ~ 1 , N Y 
RM = H 

SX = SX + (RM**2+1.) * A (M) **2 

SY = SY «■ A (H) ** 2 

IF (NV.LE.2) GO TO 712 

NVP = NY-2 

DO 708 M = 1 , N YP 

RM - M 

SX = SX - (HM + 1. ) **2*A (M) *A (M+2) 

SY - SY - A (M) * A (M + 2) 

CONTINUE 

AB2RXB = 1 2. *PI 2* (. 75*SX ♦ NU*AB2*SY) 

RYB = 1 2 . *PI 2* ( A32*3 Y + ,75*NtJ*SX) 

PHIPXX = 0, 

PHIPYY = 0. 

DO 702 M = 1 , NV 
DO 702 N = 1 , N Y 

X = B1 (M,N) *CS (M- N, XA) + B2 (M»N)*CS (M-N-2, XA) 

1 ♦ B3 (H,N) *CS (M-N+2.XA) + B4 (H , N) *CS (M+N , XA) 

2 + B5 (H ,N) *CS (M + M-2,XA) B6(N»N)*CS(M+N+2,XA) 

YMN = -16.*PI2*CS (4,YB) * X 

X = B7 (M,N) *CS (1-N, XA) ♦ B8 (M , N) *CS (H-N-2, X A) 

1 + 89 (H,N) *CS (M-M+2,XA) + BIO (M,N) *CS(fl + N»XA) 

2 + B 1 1 (M,N) *CS(M + M-2,XA) + B12 (M, N) *CS (M+N + 2,XA) 

YMN = YMN - 4.*PI2*C3 (2 , YB) *X 

RM = M 
R N - N 

X = (RM-RN) **2*3l (M,N) *CS (M-N,XA) + (Rfl-RN-2. ) **2*B2 (M, N) * 

1 CS (M-N-2, XA) * (RM-RN+2.) **2*B3 (M ,N) *CS (M-N + 2,XA) 

2 + (RM ♦ RN) * *2*B4 (M,N) *C5 (M + N , XA) + (RM + RN-2. ) **2*B5 (M, N) 

3 *CS (M + N-2, XA) + (RM+RN+2,) **2* B6 (M,N) *CS (M+N+2,XA) 

XMN = -PI 2*CS (4, YB) *X 

X = (RM-RN) **2*37 (I1,N) *CS (M-N,XA) + (RH-RN-2. ) **2*B8 (M, N) * 

1 CS (M-N-2, X A) + (RM-RN+2. ) **2*B9 (M, N) *CS (H-N+2, XA) 

2 + (RM+RN) **2*B10 (M,N) *CS (M + N,XA) + (RM + RN-2. ) **2*B 1 1 (M,N) 

3 *CS (M + N-2, XA) + (RM+RN+2. ) **2*B12 (M,N) *CS (M+N+2,X A) 

XMN = XMN - PI2*CS(2, YB) *X 

X = (RM-RN) **2*313 (M»N) *CS (M-N,XA) + (RM- RN-2. ) ** 2*B 1 4 (M , N) * 

1 CS (M-N-2, XA) + (RH-RN+2.)**2*B15(M,N)*CS(H-N+2,XA) 

2 ♦ (RM + RN) **2*B18 (M,N) *CS (M + N, XA) ♦ (RM+RN-2. ) **2*B17 (M, N) 

3 *CS (M+N-2 , XA) + (RM+RN + 2.) **2*B18(N»N)*CS(M + N+2,X A) 
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ISN 

0089 


XMN = XMN 

isn. 

0090 


PHIPXX -= 

ISN 

009 1 

702 

PHIPVY = 

ISN 

0092 


PHIPXX = 

ISN 

0093 


PHI P Y Y = 

ISN 

0094 


SIGMA X - 

ISN 

0095 


SIGMAY = 

ISN 

0096 


SIGMAX = 

ISN 

0097 

C 

SIGMAY = 
PRINT IN- 

ISN 

0098 


WRITE (6, 

ISN 

0099 


XA - XA + 

ISN 

0100 

710 

CONTINUE 

ISN 

0101 


RETURN 

ISN 

0102 


END 


- PI 2*X 

PHIPXX * XI1N*A («) *k (N) 

PHI P YY ♦ YNN* A (H) *A (N) 

12.* (1. -MU** 2) *AB2*PHIPXX 
12. * (1. -NU**2) *AB2*PHIPYY 
A32RX3 * A32*PHIPYY 
3Y3 f PrfIPXK 
3IGS1AX/12. 

SIGMAY/1 2 . 

PLANE STR EiS 33 AT POINT(XA.YB) 
701) 3IGKAX, SIGMAY 
1 . / 16 . 
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LEVEL 21.6 { MAY 72 ) 


OS/360 FORTRAN H 


ISN 0002 
ISM 0003 
ISN 0004 
ISN 0006 
ISN 0007 


COMPILER OPTIONS - NAME- H AIN , OPT=02 , LINECNT=58 , SIZE=0000K, 

SOURCE, EBCDIC, NOLIST, NODECK, LOAD, MAP, NOEDIT, ID, NOXREF 
FUNCTION 913 {A, M) 

B13 = 0. 

IF (M-N .NS* 0) 313 = -FLOAT (2*H) /FLOAT < (H-N> **3) 

RETURN 

END 
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LEVEL 21.6 ( HAY 72 ) 


OS/360 FORTRAN H 


COHPILER OPTT08S - NAME= H A t N, OPT=02, LI NECNT=58, SI2 E=OOOOK , 

SOURCE, SdCDIC, NOUST, NODECK, LOAD, HAP, NOEDIT, ID, MOXREP 
ISN 0002 PUNCTION 314 (11,11) 

ISN 0003 B 1 4 = 0. 

ISN 0004 IP (H-N-2 .NS. 0) B14 * PLOAT (H-1) /FLOAT ( (fl-N-2) **3) 

ISN 0006 RETURN 

ISN 0007 END 
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LEVEL 21.6 ( NAY 72 ) 


OS/360 FORTH AN H 


COMPILER OPTIONS - NAME- MAIN, OPT=02 , LINECNT=58,SIZE-000QIC, 

SOURCE, EiCDIC, NOLI ST, NODECK , LOAD, MAP, NOEDIT, ID, NOXEEF 
ISN 0002 FUNCTION fc15(M,N) 

ISN OOOi B1 5 « 0. 

ISN 0004 IP (M-N+-2 . NE. 0) B15 = FLOAT (fi+1 ) /FLOAT ( (M-N+2) ♦*3) 

ISN 0006 RETURN 

ISN 0007 END 
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OS/36 0 FORTH AN H 


LEVEL 21.6 ( HAY 72 ) 

COMPILER OPTIONS - N AilE= H AI K ,OPT=02 , LIN ECNT=58 ,SIZS*0000K, 

SOURCE, EBCDIC, NOLIST, NODBCK, LOAD, HAP, HOEDIT,ID, iOXREF 
I5N 0002 FUNCTION 3l6(H,N) 

ISN 0003 816 = 0. 

ISN 0004 IF (H+N ,NE. 0) 316 * FLOAT (2*M) /FLOAT ( (H+N) **3) 

ISN 0006 RETURN 

ISN 0007 END 
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OS/360 FORTRAN H 


LEVEL 21.6 ( NAT 72 ) 


ISN 3002 
ISN 0003 
ISN 0004 
ISN 0006 
ISN 0007 


COHPILER OPTIONS - NANE= N AIN, OPT=0 2, LINECNT=58, SIZ E*OOOOK, 

SOURCE, EBCDIC, NOLIST, NODECR, LOAD, HAP, NOEDIT, ID, NOIREF 
PONCTION 617 (N, N) 

B17 = 0. 

IF {H+N-2 .NB. 0) B17 - -FLOAT {H- 1 ) /FLOAT ({«♦»- 2) **3) 

RETORN 

END 
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LEVEL 21 .6 ( KAY 72 J 


OS/360 FOBTRAH fl 


ISN 0002 
ISN 0003 
ISN 0004 
ISN 0006 
ISN 0007 


COMPILER OPTIONS - NAM E= fl AIN, OPT«02 , LIHf CBT=58, SIZB*QQ00K, 

SOURCE, E8CDIC, NOLISf ,NODECK,LOAD,BAP,BOBDIT,ID,NOIBBF 
FUNCTION £ 1 8 (fl, N) 

B18 « 0. 

IF (fl*N*2 .HE. 0) B 18 « -FLOAT (K*1 ) /FLOAT ( (B+N+2) **3) 

RETURN 

END 
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